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The paper addresses the space-frequency correlations of electromagnetic waves 
in general random, bi-anisotropic media whose constitutive tensors are complex 
Hermitian matrices. The two-frequency Wigner distribution (2f-WD) for polarized 
waves is introduced to describe the space-frequency correlations and the closed form 
Wigner-Moyal equation is derived from the Maxwell equations. Two-frequency ra- 
diative transfer (2f-RT) equations are then derived from the Wigner-Moyal equation 
by using the multiple scale expansion. For the simplest isotropic medium, the result 
coincides with Chandrasekhar's transfer equation. In birefringent media, the 2f-RT 
equations take the scalar form due to the absence of depolarization. A number of 
birefringent media such as the chiral, uniaxial and gyrotropic media are examined. 
For the unpolarized wave in the isotropic medium the 2f-RT equations reduces 
to the 2f-RT equation previously derived in Part I. A similar Fokker-Planck-type 
equation is derived from the scalar 2f-RT equation for the birefringent media. © 
2008 Optical Society of America 



1. Introduction 

In Part I 6 of the series we studied the space-frequency correlation for scalar waves in random media 
as governed by the Helmholtz equation with a randomly fluctuating refractive index. To this end, 
we introduced the two-frequency Wigner distribution (2f-WD) which in the unsealed form is 



where XJ\ and U2 are the wave fields at frequencies coi and 0J2 respectively. Throughout, f denotes 
the Hermitian conjugation of vectors or matrices. The important characteristic of the definition (TjQ) 
is that the spatial argument of each wave field is scaled in proportion to the respective wavelength. 
The variables x are the so called size parameter in scattering theory when the phase velocity is 
unity. 15 
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In the weak coupling (disorder) regime we derived the two-frequency radiative transfer (2f-RT) 
equation for the two-frequency Wigner distribution. We considered several approximations, notably 
the geometrical optics and paraxial approximations. Based on the dimensional analysis of these 
asymptotic equations we obtained scaling behavior of the coherence bandwidth and coherence 
length. We also obtained the space- frequency correlation below the transport mean- free-path by 
analytically solving one of the paraxial 2f-RT equations. 

The main advantage of the 2f-RT theory over the traditional equal-time RT theory is that it 
describes not just the energetic transport but also the two space-time point mutual coherence in 
the following way. Let u(tj,Xj),j = 1,2 be the time-dependent wave field at two space-time points 
(tj,Xj),j = 1,2. Let x = (cjiXi +a;2X2)/2 and y = cjiXi — u>2^-2- Then we have 



with t = (t\ + i2)/2,T = t\ — t2,u = {uj\ + ix>2)/2,u/ = uj\ — UJ2- Here and below (•) is the 
ensemble averaging w.r.t. the medium fluctuations, * the complex conjugation and f the Hermitian 
conjugation. In comparison, the single-time correlation gives rise to the expression 



which, through spectral decomposition, determines only the central-frequency-integrated 2f WD. 
For a statistically stationary signal, ([I]) would be a function of t\ — 1-2 only. In this case different fre- 
quency components are uncorrelated and consequently only the equal-frequency WD is necessary to 
describe the two-spacetime correlation. 14 For statistically non-stationary signals the two-frequency 
cross-correlation is needed to characterize the two-spacetime correlation. 

The 2f-RT theory developed in Part I has enabled precise estimate of important physical quanti- 
ties such as the coherence length and the coherence bandwidth 6 which are medium characteristics 
relevant to communications and imaging in disordered media. 5,9 In particular, the two-frequency 
formulation is an indispensable tool for the statistical stability analysis of the time-reversal com- 
munication scheme with broadband signals in multiple-scattering media (see Ref. 5 where a 2f-RT 
equation and its solution play a key role). The 2f-RT theory developed here is expected to extend 
these results to the case of polarized waves. 

The organization of this paper is as follows. In Section [2] and Appendix A we develop the two- 
frequency formulation of the Maxwell equations for general heterogeneous dielectric in terms of 
2f-WD. In Section [3l we formulate the weak-coupling scaling limit for two- frequency Wigner-Moyal 
equation. In Section 2] we develop the multiscale expansion to find an approximate solution in 
the weak-coupling regime. In Section [5] and Appendix B, based on a solvability condition we give 
an explicit form to the 2f-RT equations for general bi-anisotropic media and in Section 15. Al we 
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derive a scalar 2f-RT equation for birefringent media. In Section I6.AI we consider the isotropic 
medium and show that the general 2f-RT equations, after a change of variable, reduces to the two- 
frequency version of Chandrasekhar's transfer equation. In Section I6.P4 16. CI and 16. D\ we examine 
three birefringent media: the chiral, the uniaxial and the gyrotropic media. In Section [7] we analyze 
the unpolarized wave in the isotropic medium in the geometrical optics regime and show that 
Chandrasekhar's equation reduces to a Fokker-Planck-type equation rigorously derivable from the 
geometrical optics of the scalar wave. 7 We derive a similar equation from the scalar 2f-RT equation 
for the birefringent media. We conclude the paper in Section [8] with a brief discussion on expressing 
the two-spacetime correlation in terms of solutions of the 2f-RT equations. 



2. Maxwell equations and Wigner-Moyal equations 

In this paper, we consider the electromagnetic wave propagation in a heterogeneous, lossless, bi- 
anisotropic dielectric medium. We assume that the scattering medium is free of charges and currents 
and start with the source-free Maxwell equations in the frequency uj domain 
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where K is, by the assumption of losslessness, a Hermitian matrix 
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with the permittivity and permeability tensors K € ,K^, and the magneto-electric tensor K x . 16 The 
Hermitian matrix K is assumed to be always invertible. 

In an isotropic dielectric, K e = el, = fil, K x = 0. In a biisotropic dielectric, K x as well 
as K e ,K M are nonzero scalars. A reciprocal chiral medium is biisotropic with purely imaginary 
K x = ix- The appearance of nonzero K x arises from the so called magnetoelectric effect. 12 Crystals 
are often naturally anisotropic, and in some media (such as liquid crystals) it is possible to induce 
anisotropy by applying e.g. an external electric field. In crystal optics, K e ,K^ are real, symmetric 
matrices and K x = 0. 2 In response to a magnetic field, some materials can have a dielectric tensor 
that is complex-Hermitian; this is called the gyrotropic effect. In general, a magnetoelectric, bi- 
anisotropic medium has a constitutive tensor ([3]) with complex Hermitian K e , and a complex 
matrix K x satisfying the Post constraint. 21 It has been shown that a moving medium, even isotropic, 
must be treated as bi-anisotropic. 4 ' 12 

In general, K is a function of the frequency uj (for dispersive media) but it turns out that if 
the frequency-dependence of K is sufficiently smooth the 2f-RT equations derived in the present 
framework have the same form as for nondispersive media; the frequency-dependence would enter 
the coefficients of the equations in the obvious way. 6 For the simplicity of presentation we shall 
assume that the medium is nondispersive. 

Writing the total field U = (E, H) we introduce the two-frequency matrix-valued Wigner dis- 
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tribution 
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where Ui and U2 are the total fields at frequencies uj\ and 002 respectively. From the definition we 
see that the variables x and p _1 have the dimension of length/time. Although the scaling factors 
in the arguments of Ui and U2 are not required for the development of the 2f-RT theory for 
the first-order (Maxwell's) equations, they are particularly useful in the case of the second-order 
(Helmholtz and paraxial wave) equations. For the consistency and continuity of presentation (see 
Section [7J) we work with the definition (JU) in the present paper. For an alternative development 
of the 2f-RT theory for Maxwell's equations in terms of the 2f-WD without the scaling factors, we 
refer the reader to Ref. 8 

First note the symmetry of the Wigner distribution matrix 



W f (x,p;u;i,u;2) = W(x, p; lo 2 , wi). 



(5) 



In other words, the right hand side of @ is invariant under the simultaneous transformations of 
Hermitian conjugation f and frequency exchange lo\ <-> o>2. 

In what follows we shall omit writing the arguments of any fields if there is no risk of confusion. 

We put eq. ([2]) in the form of general symmetric hyperbolic system 17 
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The matrices iTj,j = 1,2,3 are related to the photon spin matrices. 1 

Throughout this paper the dot notation, is used exclusively for directional derivative as 
in p • V = Pjd Xj . All vectors are treated as matrices and the scalar product is just the matrix 
multiplication between row and column vectors. All vectors are taken to be, by default, column 
vectors, unless explicitly transposed. 

Applying the operator Rjd/dxj to W and using flS} we obtain 
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whose derivation is given in Appendix A. From ([7]) and ((5|) we also have 



2iWpjRj -2i j W(x,p + 



(8) 



2lo 2 
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dq. 



Here and below K stands for the Fourier transform (spectral density) of K as in 

K(x) = J e ixtq K(q)dq. 
For a Hermitian K we have K(p) = K'(— p), Vp. 
3. Weak-coupling limit 

As in Part I 6 we consider the weak coupling regime with the tensor 

K(x) = K (l + V / lv(j)), £«1 (9) 

where the Hermitian matrix Kq represents the uniform background medium and y/iV represents the 
relative fluctuations of the permittivity-permeability tensor. The small parameter £ describes the 
ratio of the scale of the medium fluctuation to the propagation distance. In an isotropic dielectric, 



Kn 



where e and p, are electric and magnetic susceptibility, respectively. In general Ko is a Hermitian 
matrix and its blocks, as in ([3|), are denoted by Kq, Kq , Kq , Kq^, respectively. To preserve the 
Hermicity of K and Kq the matrix V must satisfy 
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We shall assume below that Ko is either positive or negative definite. Otherwise, the materials 
would be lossy since the refractive index is not real-valued if Ko is not sign-definite. A negative- 
definite Ko gives rise to negative refractive index which is a hot topic in metamaterial research. 18-20 
To fix the idea, let us take Ko to be positive definite. With minor notational change, our method 
applies equally well to the negative definite case. 

We assume that V = [Vij] is a statistically homogeneous random field with the spectral density 
tensors * = [$y mr J, * = [*i? mn ] such that 



(^(x^Jy)) 
This implies the following relations 
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In the case of real- valued V, $ = 1 J r . The spectral density tensors have the basic symmetry 



^mn(p) = *mnij(p), ( 15 ) 
^ijmn(-p) = *mnii(p)> (16) 



Furthermore, eq. (jlOp implies that 



Ko,ij^mnjl(p) = Ko,lj®mnji{p) (17) 
K ,ij^mnjl(p) = KZyVmnjifa) (18) 

As in Part I, we consider the regime where the wavelengths are of the same order of magnitude as 
the correlation length of the medium fluctuations by rescaling the frequencies u>« = u>j/£,j = 1,2. 
This choice of frequency scaling results in strong scattering by the medium heterogeneities. For 
ease of notation, we drop the tilde in ojj below. To capture the high frequency behavior of the wave 
field we redefine the 2f-WD as 

W(x,p) (19) 

W e -ip< yUl (*+i>>t(^-^)<(y. 
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We also assume that ui,u)2 — > u as £ — > such that 
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with a fixed constant (3. The governing equations for fj 191) become 

R ^ W = -TPi R ;W + ^K W (21) 
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where x = x/-£ is the fast spatial variable. In order to cancel the background effect we multiply eq. 
Plj) by Kq" 1 from left, ([22]) by K^ 1 fr om right and add them to obtain the symmetrical form 

K ° lR ^ W + ^- WRjK °" 1 (23) 
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This is the equation that we shall work with to derive the 2f-RT equations emplying the multiscale 
expansion (MSE). 6 ' 17 Note that eq. (|23p is invariant under the simultaneous transformations of 
Hermitian conjugation f and frequency exchange uj\ <-> a>2- 

If, instead of adding the two equations, we subtract them then we obtain the anti-symmetric 
form 

"T W + K o 'B^W - ^-WR,K 1 



dxj dxj 



+ T [K - 1 p,R J W + Wp J R j K 
21 



o 1 ] 



2% 

- 7f/[ e,,,i/ "^ (q)w(p -^ ) 

+W(p- -5L)vt(q) e - iqti/<JJ2 dq. (24) 

2tJ2 J 

Eq. (|24p requires a different treatment and will not be pursued here. However, the leading order 
l~ x terms of eq. (|24p impose a constraint which will be discussed in the Conclusion. 

4. Multiscale expansion 

The key point of MSE is to separate the fast variable x from the slow variable x and make the 
substitution 

R.cL.W -> R— W + ^Rj— W 

ct.WR, -> — WRj+r 1 — WR r 

The idea is that for sufficiently small £ the two widely separated scales, represented by x and x 
respectively, become mathematically (but not physically) independent. 

We posit the expansion W = W + \/fWi + W2 + substitute it into eq. (|23|) and equate 
terms of same order of magnitude. 

4- A. Leading term 
The £ _1 -terms yield 

K ° lR ^ W + J- WRjK ° 1 (25) 
+2i [K-^-RjW - WpjRjK^ 1 ] = 0. 

We hypothesize that the leading order term W = W(x, p) be independent of the fast variable x. 
Thus the first two terms of (125j) vanish so the equation reduces to 



K^ Pj KjW - W Pj RjK^ = 0. (26) 

Eq. (|26p arises also in the equal-time RT theory 17 and can be solved as follows. For a positive (or 
negative) definite Kq, consider the eigenvalues and eigenvectors {e CT,a } of the matrix Kq 1 pjHj 
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where the index a keeps track of the multiplicity and hence depends on a. As Kq~ PjR,- is Hermitian 
with respect to the scalar product defined by a^Kob,Va, b G C 6 , the eigenvalues are real and the 
eigenvectors form a complete set of Ko-orthogonal basis in C 6 . Alternatively, we may work with the 

— 1/2 — 1/2 

Hermitian matrix K Q PjRjKq in the image space, with the standard scalar product, under the 
transformation Kq^ 2 . Let the eigenvectors {e°"' a } be normalized such that e <T,a ^'Koe T '^ = S a>T 6 a ^. 

Clearly, the eigenvalues f2 CT as a function of the wavevector p define the dispersion relations. 
For general bianisotropic dielectric, it is easy to check that Cl° = is always an eigenvalue with 
eigenvectors 

e^Ctf-Q, e°. 2 (p)~Q. (27) 

Since Ko is invertible, it follows that the null space of Kg 1 pjHj is spanned by these two non- 
propagating modes. It is easy to check that {d°"' Q t(p) : d°" ,a (p) = Koe°"' Q (p)} are the left eigen- 
vectors of K^pjKj and {d ff ' Q (p)}, {e r ^(p)} are co-orthogonal with respect to the standard scalar 
product: 

d*.«t(p)e r '<(p) = $ ff , T <W- ( 28 ) 

This relation will be useful in deriving the 2f-RT equations (see Appendix B). 

Throughout the English indices represent the spatial degrees of freedom while the Greek indices 
represent the modal and polarization degrees of freedom. It is important to keep this distinction in 
mind in the subsequent analysis. The Einstein summation convention and the Hermitian conjuga- 
tion are used only on the English indices. 

It can be checked easily that the general solution to (|26j) is given by 



W(x,p) = £W£(x,p)E ff '°<(p,p) (29) 

a,a,C 

where are generally complex-valued functions and 

E^p, q) = e ' a (p)e ' c \q). (30) 

Likewise we define 

D^(p,q)=d^(p)d^ t (q). (31) 

The linear span of {E T,Q ^(p, p), Vr, a, £, p} is a Hilbert space, denoted by 9JT P , for each p / 
with the scalar product Tr [H+KGK] , H, G £ 9K p . The matrices W CT = [W^], free of the English 
indices, are called the coherence matrices. 

For x-independent W, the constraint that the electric displacement and the magnetic induction 
are both divergence-free yields, on the macroscopic scale, 

(±V,±V) K W = 
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which, in view of the definition (|19p . is equivalent to 

(±p t ,±p t )K W(x,p) =0. 



(32) 



Hence by ([2~T|) d ' 5 'W = and by ([28]) W° = where W° in (f29|) is the coherence matrix associated 
with the non-propagating mode f2° = 0. 

4-B. Correctors 

The £ -1 / 2 -terms yields the equation 



+2i [Ko^jRjW! - WiPjRjKq x ] 

= 2i [ dq. e iqt * M V(q)W(p - -5-) 
J [ 2u\ 



-W(p - -5-)Vt(q) e - 4qt ^ 2 
2w 2 



(33) 



where, as in Part I, 6 we have added a small regularization term. The reader is referred to Part I 6 for 
the discussion of the choice of the regularization parameter. Physically, the sign of the parameter 
(positive here) amounts to choosing the direction of causality. 
We Fourier transform eq. (|33|) in x 



and posit the solution 



-*2fWi(k,p) + Ko 1 %R i Wi(k, p) +Wi(k,p)/c,-R,-Ko 1 
+2 [KoVjRjWx^p) - W 1 (k,p) Pj K j K 1 

2 [v(wik)W(p - ~) - W(p + |)vt(-u; 2 k) 



Wi(k,p) = £ Q(k,p)E^(p + |,p- |) 



(34) 



(35) 



where are generally complex numbers. Note that the two arguments of E°" ,Q ^ in (|35|) are at 
different momenta p + k/2,p — k/2. 

We substitute and ((35]) into eq. {33]) and multiply it with d <T '° t (p+ 1) from the left and with 
d CT '^ (p — ^ ) from the right and solve the resulting equation algebraically. This yields the coefficients 

C£ c (k,p) 

jr(p + k)-fr( P -k)- 
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d ^t (p+ )v( Wlk )W CT c (p 
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When the leading term W is invariant under the simultaneous transformations of Hermitian con- 
jugation f and frequency exchange lj± <-> so is Wi. This invariance is manifest in the relation 

C<a(-k,P;^i,^2) = C^ c (k,p;cj 2 ,wi)- 
Finally the 0(l)-terms yield the equation after adding a regularizing term 2W2 

2£W 2 + K^K~W 2 + ^W 2 R i K^ 1 

+2i [Kq 1 pjB. i W 2 - Wapj-Rj-Ko x ] = F (37) 

with 

F = 2i J 



eiq txM v(q)Wi(p _^_ ; 

lull 



-W 1 (p--5-)Vt(q)e^ ti M 

2oJ 2 

K,: ; H,^-W - ^WR^ 1 . (38) 

Note again that F is invariant under the simultaneous transformations of Hermitian conjugation 
f and frequency exchange uj\ <-> uj 2 . We can, but need not, solve eq. ([37]) explicitly as eq. (|3"4"]l . 
However, in order for the second perturbation -fW 2 to vanish in the limit t — > 0, F must satisfy the 
solvability condition 

Jim Tr ^G^KoFKo^) = (39) 



for all random stationary matrices G satisfying eq. (|25p . This can be seen by transforming eq. ([3] 
into Tr (G^KodSTjKo) which by eq. (ESJ) implies 2£Tr (Gl~KoW 2 Ko) = Tr <GtK FK > and hence 

432)- 

Fortunately, we do not need to work with the full solvability condition (|39p . It suffices to demand 
(j39|) to be fulfilled by all deterministic G, independent of x, such that 

K^pjKjG - GpjRjKq 1 = 0. (40) 

In other words, as in ([29]) . we consider only a subspace of the solution space of eq. (f25j) and replace 
(39D by 

lim Tr (D^ f (p, p) (F(x, x, P ))) = 0, Vr, £, u, p, x, x (41) 

where D T '^ are defined in (|3~TI) . As noted above, ([23]) . (j33j) and (f38l) are invariant under the 
simultaneous transformations of Hermitian conjugation f and frequency exchange u)\ <-> u; 2 and 
therefore eq. (|4ip must also be invariant under the same transformations. 

To summarize, we have constructed the three-term expansion W + v^Wi + £W 2 which is an 
approximate solution of the 2f Wigner-Moyal equation in the sense the left hand side of (|23p 
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subtracted by the right hand side of (|23|) equals exactly 

VI[ - 2Wi + K ; R ; ^ W ; + i),. W.KjK,, 1 
-2iV2 [ [e 4qti /^V(q)W 2 (p - -5-) 

dq 



1 \-Ti-tf„ , \ -Mi t x/w2 



- W2 ( p -2^) VT(q)g 



+£ 



2W 2 + Ko ^jflx.Wa + d Xi W 2 R,Ko 1 



which vanishes in a suitable sense as t —> 0. 6 

With (|35]) - (f36|) and (f38l) . eq. (jITj) is an implicit form of the 2f-RT equations that determines the 
leading order coherence matrix. Our next step is to write (|41j) explicitly in terms of the spectral 
densities of the medium fluctuations. 

5. 2f-RT equations 

Calculation with the left hand side of eq. (|4ip is tedious but straightforward as it involves only the 
second order correlations of V. This is carried out in Appendix B. 

To state the full result in a concise form, let us introduce the following quantities. Define the 
scattering tensors © r (p, q) = [5^ a ^(p,q)] as 



S£«<(p,q) 

= <c*(pK'° (q) W(p - q))4''(p)^' c *(q) 

Using (|15|) -(fT8j) one can derive the alternative expressions for S: 



(42) 



= e^*(p)d^(q)** /ss ; ( w (q - p))<^(p) e ^*(q) 

= ^(p)er(q)$ /9S1 ;( W (q-p))e^(p)rff(q). (43) 
With (fl5|) . (fTBT) and (l43j) it is also straightforward to check that 

S£««(P,q) =^ aK (p,q) =5 f ^(q,p) (44) 
For any 9Jl p -valued field G(p) define the (£, ^-component of the tensor 6 T (p, q) : G(q) as 

[6 T ( P) q) : G(q)]^ = ^<S^ c (p,q)G aC (q). 

Define the tensors 5] r = J analogous to the total scattering cross section as 

s t (p) = nj 5(V( P ) - n T (q))© T (p,q) = Mq (45) 

-i^ (^(p)-^(q))- 1 6 T (p,q) :Idq. 
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The 2f-RT equation then reads as 

V p ft T • V X W T = 2vrw 3 J s(p, T (p) - O r (q)) (46) 

xe- i/3(q - p)tx S T (p,q) : W T (q)dq 
-u, 3 (V(p)W T (p) + W T (p)5Tt( p ) 



Vr. 



Introducing the new quantity 



2H T = e- i/3ptx W T (p) 
we recast eq. (|46p into the following form 

V P U T • V x 2H r + i/Sp • V p fi T 2B T (47) 
= 2itlu 3 J <5(V(p) - fi T (q))e(p,q) : 2H r (q)fiq 

s-( P )2ir(p) + 2ir(p)s-t(p) 



-a; 3 



This is the Rayleigh-type scaling behavior typical of a random dielectric. The cubic, instead of 
quartic, power in uj is due to the appearance of uj as the scaling factor in the definition of 2f-WD 
(JT9J) . The quartic-in-w law is recovered upon replacing x by x/w on the left hand side of (|47|) . 



5. A. Decoupling: scalar 2f-RT equation 

Although, in view of (I27j) . the zero eigenvalue Q° = has multiplicity two in general, the nonzero 
eigenvalues in media other than the simplest isotropic medium often have multiplicity one as we 
shall see in Section [6l This is closely related to the birefringence effect. Under such circumstances, 
the 2f-RT equations take a much simplified form which we now state. 

Because j = 1,2,3,4 are simple (multiplicity one), expression (|29l) reduces to 

W(x,p) = ]T^(x,p)E CT (p,p). 

a 

In other words, the coherence matrices become scalars and the different polarization modes decou- 
ple. Consequently (|4"6j) becomes a scalar equation 

VpfT • V X W T = 2ttuj 3 j S(n r (p) - ft T (q)) (48) 

xe -i/3(q-p)tx 6 r q) ^r ( q)dq 

-2lu 3 Z t ( P )W t (p), Vt 

where 

6 r (p,q) = C(p)e[(q)^/ 5 (^(p-q))dKp) e r(q) ( 49 ) 
E T (p) = ^yj(^(p)-^(q))6 T (p,q)dq. (50) 

Note that the Cauchy singular integral term in (|50p disappears whenever S T and W commute as in 
the scalar case. Prom (|48p we can derive the scalar equation for the quantity 2B r = e~ t ^ p '' x W T (p) 
as before. 
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6. Special media 

In this section, we consider the eigenstructure of the dispersion matrix Kq 1 pj^-j associated with 
the various background media for which the scattering tensor can be computed explicitly. 

6. A. Isotropic medium 

In the simplest case of an isotropic medium, there are two nonzero eigenvalues: il + (p) = 
co|p|,fl~(p) = — Co | p | , each of multiplicity two. Let p = p/|p| and let p^pj be any pair of 
unit vectors orthogonal to each other and to p so that {pi", pj, p} form a right-handed coordinate 
frame. Let {q"t,q7,q} be similarly defined. The eigenvectors are 



e+'+(p) 




(P) 




■(P) 




IP J 




Denote the spectral densities of e and p, by <l? e and respectively and denote the cross spectral 
densities by <J? e/i , (i?^. We have & T = [S^ ua ^\ with 



«S£*<(p,q) 



d> e (u;(p - ^tfMl^pl 
-fc^Mp - q))pi t qlq: Ct p:" 

+^( W (p- q ))p-« t q I - q - Ct p I - 



for r, £, a, £, v = ±. Eq. (|46p can now be written as 



c p • V x W d 



7ru; 3 |p| 2 



2 / e- i/3 ( q - p ) tx 5(|p 



4c 

x6(p,q) : W ± (q)dq 

- |q|)s(p,q) : MqW 



|q| 



(51) 



(52) 



± 



<m IpI - |q| )©(p,q) : Id q 



-w= 



The property ([4"4"]1 and the expression (f5Tj) imply that SJ avi -{p, q) = S^ ai/ ^(q, p) and hence the 
Cauchy principal value term in (|45|) disappears. 

Often, in a scattering atmosphere for instance, jl = is a good approximation and in such case 
the only nonzero term in the scattering kernel is 

1 



5 [m<(p>q) 



-^Hp-q^pjqiq^pl 



(53) 



This is the setting for which S. Chandrasekhar originally derived his famous equation of transfer 3 
and eq. (|52p is just the two-frequency version of Chandrasekhar's equation. 

In the same setting, the new features in (|52p beyond Chandrasekhar's transfer equation are 
the frequency shift (3 and the general form of the power spectrum $ e . In Chandrasekhar's and 
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other cases, the medium consists of randomly distributed particles of smaller size than the wave- 
length. 10,15 Such a discrete medium correspond to a random field V that is a sum of <5-like functions 
randomly distributed according to the Poisson point process whose spectral density tensor 3> can 
be calculated. 



6.B. Chiral media 

A chiral medium is a reciprocal, biisotropic medium with the constitutive matrix 



e I ixl 



where x 6 ^ is the magneto-electric coefficient. To maintain a positive-definite Ko we assume 
X 2 < eofJ-o- We then have 



K l p j 'R j 



co 



zl —ikI 




-px 


inl z _1 I 




px 



(54) 



where z = \J '//qAo > is the impedance and k = x c o is the chirality parameter. The four non- 



zero simple eigenvalues are fi 1 = co|p|(l + k) 1 ,f2 2 = co|p|(l 
c o I P I ( — K ~ I) -1 an d their corresponding eigenvectors are 



= cq\p\(k- iy\n 4 



6 ~ _pi_,pi 



e 2 ~ 



2 Z / 



e 3 ~ 



-«pi + pi 
£l + i pi 

Z 2 / 



4 

, e ~ 



^Pi + p! n 

\ Z 1 2 / 



Note also that Q 4 = — fi 1 ,^ 3 = — f2 2 . As |«| < 1 (since x 2 < ^0/^0)1 e 1 ^ 2 are the forward propa- 
gating modes and e 3 ,e 4 the backward propagating modes. 

For the medium fluctuation V we may use the biisotropy form 

a\ i/iobl 
-ieobl al 

where a,b £M are stationary random functions of x with power spectral densities & a ,&t> and the 
cross-spectral density $ a h. This particular form is derived from the commutativity relation (|10p . 

The splitting into two distinct positive dispersion relations is a case of birefrigence where two 
distinct phase velocities, co/(l ± k), arise depending on the polarization. As discussed in Section 
l5.Al due to the birefringence the chiral medium does not depolarize the electromagnetic waves. For 
the sake of space, we leave to the reader to work out the scattering tensor from (|49 |) -(|50 [) . 

6.C. Birefrigence in anisotropic crystals 

Generally speaking, an anisotropic medium permits two monochromatic plane waves with two 
different linear polarizations and two different velocities to propagate in any given direction. 2 This 
again gives rise to the birefringence effect. 

The only optically isotropic crystal is the cubic crystal. In the system of principal dielectric axes, 
the permitivity-permeability tensor of a crystal, which is always a real, symmetric matrix, can be 
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diagonalized as Ko = diag[e x , e y , e z , 1, 1, 1]. One type of anisotropic crystals are the uniaxial crystals 
for which e x = e y = e± ^ e z = en (if the distinguished direction, the optic axis, is taken as the 
2- axis). There exist two distinct dispersion relations for the forward modes 



n° =JjeL, n« =M + pI±A, 



The backward modes correspond to —$7°, — f} e . The corresponding wavevector surface consists of a 
sphere and an ovaloid, a surface of revolution. £1° corresponds to the ordinary wave with a velocity 
independent of the wavevector while £l e corresponds to the extraordinary wave with a velocity 
depending on the angle between the wavevector and the optic axis. 2 

Let d°,d e be the associated left eigenvectors. Set K = diag[ej_, ej_, en] and let a a solve the 
following symmetric eigenvalue problem: 

- p x (Kg) -1 p x a CT = (1T)V, a = e,o. (55) 

Then the left eigenvectors d CT can be written as 

The same formula applies to the backward modes. Eq. (|55p has the following solutions 

2 2 

a e = (-p 2 ,Pi,0)t, a ° = (p 1 , P2 - P l±Pl)\ 

P3 

from which we see that the wave is linearly polarized. 

The other type of anisotropic crystals is the biaxial crystals for which there are also two distinct, 
but more complicated, dispersion relations, both associated with the extraordinary waves. 2 In 
contrast, the two distinct dispersion relations of a chiral medium give rise to two ordinary waves 
as the two wavevector surface consists of two concentric spheres centered at p = 0. 

It should be emphasized that a plane wave propagating in an anisotropic crystal is linearly 
polarized in certain planes whereas a plane wave propagating in the isotropic medium is in general 
elliptically polarized, and is linearly polarized only in particular cases. In the anisotropic as well 
as the chiral media, the different polarizations decouple in the radiative transfer equations and the 
depolarization effect is absent. 

6.D. Gyrotropic media 

In the presence of a static external magnetic field H cxt the permittivity tensor Kq is no longer sym- 
metrical; it is generally a complex Hermitian matrix. Here we consider the simplest such constitutive 
relation 

D = e E - ig x E, B = H (57) 
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where g = /H cxt , / G K, is the gyration vector. Equivalently, we can write 

E = 2 1 i9 ( £ oD + tg x D - -gg f D ) . 

4 ~ |g| 2 V e / 

In this case there are two distinct forward dispersion relations 12 

n =cq|p + — g|, o =c |p-— g| 

where cq = 1/ yfe~Q. Clearly the wave- vector surface consists of two spheres of the same radius but 
different centers. This should be contrasted with the case of chiral media for which the wave- vector 
surface consists of two concentric spheres of different radii. 

The associated (left) eigenvectors d a , a = 1,2 can be written as in (|56j) with a CT solving (|55f) and 
Kq corresponding to (|57|) . Let g = g\p]_ + g2P± + 93P- We can write the three-dimensional vector 
a°" as a°" = + 'jap'x with 



,. _ g 2 2 -gi 2 -(-i)V(^ + g 2 2 ) 2 + 4^ „_ 12 

2(5i92 - «e 53) 

We see that the wave is in general elliptically polarized or linearly polarized when g is orthogonal 
to the wavevector p and circularly polarized when g is parallel to p. Again, the simplicity of the 
eigenvalues implies that depolarization is absent in the gyrotropic media. 

7. Geometrical 2f-RT 

We have seen in Section T5.AI how a scalar 2f-RT equation naturally arises in a birefringent medium. 
In this section, we show that a scalar 2f-RT equation can also arise in a depolarizing medium such as 
the isotropic medium discussed in Section fc.Ai Depolarization can mix different polarization modes 
and result in scalar-like coherence matrices W r ~ W t I,t = ± (see Section [6. Al for notation). The 
other purpose of this Section is to show smooth transition from (|46|) to the Fokker-Planck equation, 
previously derived for the scalar waves, 6 ' 7 in the geometrical optics through rapid depolarization. 

Let us start with the general setting and replace <&(•) in (|47|) by 7~ 4 $(-/7) where the small 
parameter 7 is roughly the ratio of the wavelength to the correlation length of the medium fluctua- 
tions. In other words, we consider the geometrical optics regime. The quartic power in 7 is indicative 
of the Rayleigh-type scattering. Consider the change of variable q = p + 7k in the scattering term 
of (|47p . With this and the ansatz 2H r = e _ *^ p t x W r the scattering term becomes approximately 

2WV 1 J dk5(V(p + 7k) -ft T (p))c^*(p)^(p) 
x $ si/g ( W k) + 7k) e T g a * (P + 7k) 

a 

x 2U T (p) + 7 k • V p 2B T (p) + -Lkikjd^fnTfp)] . (58) 
The first term in (|58|) cancels exactly with 5] T (p)2n r (p) + W T (p)'S T ^(p) on the right hand side of 
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TTOJ 



([37|1 . The second term in ([58]) yields the first order differential operator 

■^d^^^^^W (59) 

a 

x y fcj&j5'(k- V p ft T )$ si / fl (c<;k)kdk 

a 

krfik ■ V p nA$ sifg (ujk)kdk 



• v p *m r 

where 5' is the derivative of the Dirac-delta function. And the third term in (|58p yields the second 
order differential operator 

7ru 3 d T /*d T f ' v J2 e i ae 9 a * [ <*(k-V p fi T ) 

x<Z> sifg (ujk)k m k n dkd Pm d Pn W T . (60) 

In order to match the left hand side of (|47p which is a scalar in case of complete depolarization, 
(|60p and each term in (|59p must be proportional to as well. This happens, for instance, for the 
isotropic medium with ([53]) . In this case, 



flSI = ttc^^- / k- V p (k-p)<5'(k-p)$ £ (a;k)kdk- V p 2B T , 

3 

m = 7t5^ u ^J s(k-p^ e (uk)k m k n dkd Pm d Pn mr 

and hence the 2f-RT equation (pT7|) becomes 

±c p- ± z/3c |p|2n ± = Vp-DVp^ 

with the diffusion coefficient 

D = J 5[k ■ pj$ e (a;k)kkt£ik (61) 

which is the same Fokker-Planck-type equation derived by a rigorous, probabilistic method from 
the geometrical optics of the scalar wave previously. 

Applying the same procedure to the scalar 2f-RT equation for the quantity 2U T of the birefringence 
case discussed in Section I5.A| we obtain 

VpJT • V X 2H T + ipp • V p fT2ir = V p • DV P 2B T 

where the diffusion coefficient D is given by 



D( P ) = ™ 3 C(pK(pK(pK*(p) 



J 5(k--p)$ sifg (uk)kkUk. 
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8. Conclusion 

Starting with the symmetrical Wigner-Moyal equation (|23|) . we have systematically derived the 
2f-RT equations (|46p . (|47p and f|48j) in the radiative transfer regime characterized by the weak 
coupling scaling ([9]) . The main assumptions on the medium are that the background is uniform and 
has a either positive or negative definite constitutive matrix and that the fluctuations are zero-mean 
statistically homogeneous processes. 

We not turn to the antisymmetrical Wigner-Moyal equation (|24[) and discuss the consequence of 
its leading order terms which are 



In view of ()29[) this is equivalent to 1 = f2 CT (p). Note again the variable p has the dimension of 
inverse velocity. Therefore the two-spacetime correlations of the time dependent polarized wave 
field u are given approximately by 



with x = cj(xi + X2)/2,y = w(xi — X2)/£ where 2U 7 " = [2Jjy is the solution to eq. ([4"7]) and dfi(p) 
is the area element of the surface fi CT (p) = 1. 

Parallel to the case of scalar waves one can also work out the implications of the polarization on 
the problems of imaging and time-reversal communications, as discussed in the Introduction and 
references therein, from the 2f-RT theory developed here. 

9. Acknowledgements 

Research is supported by ONR Grant N00014-02-1-0090, Darpa Grant N00014-02-1-0603. 
References 

1. I. Bialynicki-Birula, "Photon wave function," Prog. Opt. 36, 245-294 (1996). 

2. M. Born and W. Wolf, Principles of Optics, 7-th (expanded) edition (Cambridge University 
Press, 1999). 

3. S. Chandrasekhar, Radiative Transfer (Dover Publications, New York, 1960). 

4. D. K. Cheng and 1. A. Kong, "Covariant descriptions of bianisotropic media," Proc. IEEE 56, 
pp. 248-251 (1968). 

5. A. C. Fannjiang, "Information Transfer in Disordered Media by Broadband Time Reversal: 
Stability, Resolution and Capacity," Nonlinearity 19 (2006) 2425-2439. 

6. A.C. Fannjiang, "Two-frequency radiative transfer and asymptotic solution," J. Opt. Soc. Am. 
24 (2007), 2248-2256. 

7. A. C. Fannjiang, "Space- frequency correlation of classical waves in disordered media: high- 
frequency and small scale asymptotics, " Europhys. Lett. 80 (2007), 14005. 



2W = K^pjKjW + Wpj-RfK, 



— l 
o 




(62) 



18 



8. A. C. Fannjiang, " Mutual coherence of polarized light in disordered media: Two-frequency 
formulation extended", To appear J. Phys. A: Math. Theo., arXiv:0708.4044, 

9. A.C. Fannjiang and P.M. Yan, "Multi-frequency imaging of multiple targets in Rician fading 
channels: stability and resolution," Inverse Problems 23 (2007) 1801-1819. 

10. A. Kokhanovsky, Polarization Optics of Random Media (Springer, 2003). 

11. H.N. Kritikos, N. Engheta and D. L. Jaggard, "Symmetry in electromagnetic media: A spinor 
viewpoint," in Electromagnetic symmetry C.E. Baum and H. N. Kritikos, eds. (Washington, 
D.C. : Taylor & Francis, 1995), pp 185-230. 

12. L.D. Landau, E.M. Lifshitz and L.P. Pitaevskii, Electrodynamics of Continuous Media (Elsevier 
Butterworth-Heinemann, Oxford, 1984) 

13. Lindell I V, Sihvola A H, Tretyakov S A and Viitanen A J, Electromagnetic Waves in Chiral 
and Bi-isotropic Media (Boston, MA: Artech House, 1994) 

14. L. Mandel and E. Wolf, Optical Coherence and Quantum Optics (Cambridge University Press, 
1995). 

15. M. Mishchenko, L. Travis, A. Lacis, Multiple Scattering of Light by Particles: Radiative Transfer 
and Coherent Backscattering (Cambridge University Press, Cambridge, 2006). 

16. T.H. O'Dell, The Electrodynamics of Magneto-Electric Media (Amsterdam: North-Holland, 
1970). 

17. L. Ryzhik, G. Papanicolaou and J.B. Keller, "Transport equations for elastic and other waves 
in random media, " Wave Motion 24, 327-370 (1996). 

18. D.R. Smith, W. Padilla, D.C. Vier, S.C. Nemat-Nasser and S. Schultz, "A composite medium 
with simultaneously negative permeability and permittivity," Phys. Rev. Lett. 84, 4184-4187 
(2000). 

19. D.R. Smith, J.B. Pendry and M.C.K. Wiltshire, " Metamaterials and negative refractive index, 
" Science 305, 788-792 (2004). 

20. V. G. Veselago, "The electrodynamics of substances with simultaneously negative values of e 
and /i," Sov. Phys. Usp.10, 509 -514 (1968). 

21. W.S. Weiglhofer and A. Lakhtakia, Introduction to Complex Mediums for Optics and Electro- 
magnetics. Bellingham, WA: SPIE Press, 2003. 



19 



A. Derivation of Wigner-Moyal equation 

Applying the operator Rjd/dxj to W we have 

R^wfx.p) 



/^ P ' yR 4 U < + - & dy 



(2vr) 3 y J V 2oo 1 / dyj ZK u 2 2u; 2 " 

Integrating by parts with the second integral and using ([6]) we obtain 

R^W(x,p) 

_^L p R f e -P t y Ul (- + JL)uJ(— - -^W. 

(2tt) 3 ^ J 7 U wi 2cV 2V u 2 2<V * 

Inserting the spectral representation of K intro the equation and using the definition (|19|) we then 
obtain ©. 

B. Calculation of eq. (14 1|) 

B.A. Propagation terms 
We first show that 



Tr 



D r ^t( K - 1 R J ^W + ^WRiKo 1 ] 



2v p n T • v x w^. 



Consider the following calculation 



K„ ; Rji),. W 

VptKo^Ry] • [V x W%]e^e^t 
V p [K^ Pj R je ^] • [V x W£]e<*t 
-K-^-R,- [V p e^] • [V x W£]e^ 
V p fT • V x V^ c e CT ' a e CT '^ 
+ (n° - K -^R,) [V p e*'«] • [VxWScK'tf 



Upon the operation Tr[D T, ^( • )1 the second term vanishes while the first term reduces to V p f2 T ■ 
V X WX, by (|28p and the fact that d T ^ is a left eigenvector of the matrix Kq Pj'Rj with the eigenvalue 

sr. 

The other term Tr [D r '^ d Xj WRj Kq 1 ] gives the identical result. 
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B.B. Scattering kernel 

The (s, j ) -element of the matrix 

( f > +k » r 4 w+ ^ ¥ » 

= 2i J dq/e^^W^-^-) 
-W 1 (p-^-)Vt(q) e ^ t -/-^ 

has the expression 

2*4 J dk(^ CT ( P + k) -o ct ( p ) -ay 1 

xd}> a *(p + k)$ /9si ( Wl k)f; ( (p)e^(p)^(p + k,p) 

-2*4 [ dk (V(p + ^(— + l)k) - CT (p + - l)k) - u 

J \ 2 L0\ 2 L0\ 

xe i(1 -- )kti t;(p + 5(1 + £)k)*} ff-i (-u*k) 

xer(p+ i (1 + ^ )k) ^ (p+ ' ( ^_ 1)k) 
xE , T<(p+ > (1 + ^ )k , P+ ^_ 1)k) 

-2»wl / dk (V( P + — )k) - o CT ( P - + i)k) - 

xe* (1 -^ )kt ^ ? (P - + l)k)$ /9jn ( Wl k) 

x^*(p + 5(1 - ^)k)e™(p - 1(^1 + l)k) 
x^T C (P + ^(l--)k,P-^ + l)k) 

+2^ 1 dk (ir (p) - n-( P - k) - rf)" 1 w^( P )^ fgjn (-u 2 k) 

xe^*(p)^(p-k)£;r C (P,P-k). 
Using the identity 



the symmetry properties ([T^]) - (fTU]) and (|2T|) we obtain in the limit £ — ► eq. from 
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